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ABSTRACT 

This paper presents the results of collisional evolution calculations for the Kuiper 
belt starting from an initial size distribution similar to that produced by accretion 
simulations of that region - a steep power-law large object size distribution that breaks 
to a shallower slope at r ~ 1 — 2 km, with collisional equilibrium achieved for objects 
r < 0.5 km. We find that the break from the steep large object power-law causes a 
divot, or depletion of objects at r ~ 10 — 20 km, which in-turn greatly reduces the 
disruption rate of objects with r > 25 — 50 km, preserving the steep power-law behavior 
for objects at this size. Our calculations demonstrate that the roll-over observed in the 
Kuiper belt size distribution is naturally explained as an edge of a divot in the size 
distribution; the radius at which the size distribution transitions away from the power- 
law, and the shape of the divot from our simulations are consistent with the size of the 
observed roll-over, and size distribution for smaller bodies. Both the kink radius and 
the radius of the divot center depend on the strength scaling law in the gravity regime 
for Kuiper belt objects. These simulations suggest that the sky density of r ~ 1 km 
objects is ~ 10 6 — 10 7 objects per square degree. A detection of the divot in the size 
distribution would provide a measure of the strength of large Kuiper belt objects, and 
constrain the shape of the size distribution at the end of accretion in the Kuiper belt. 



1. Introduction 



The Kuiper belt is a population of planetesimals outsi de the orbit of Neptune which ex hibits 



high inclinations and eccentricities by many belt members (ITruiillo et a" 



has a very low mass, M < 0.1 M ffi (|Gladman et al 



2001 



Bernstein et al 



20081 ). The high relative encounte r velocities, v re i ~ 1 km s 1 (Dell'Oro et al 



collisions of the largest members (|Davis &: Farinella 



1997 



2001 



2004; 



Durda k, Stern 



Brown 



2001 



and 



Fuentes fe Holman 



2001 



200' 



) and infrequent 
make the growth 

of Eris-sized bodies impossible over the age of the Solar system. Accretion in the early stages of 
planet-b uilding must have been in a more dense and quiescent environment allowing large objects 
to grow dstern & Colwelllll997h . 



The size distribution of the Kuiper belt is the result of the accretionary and collisional processes 
that have gone on in that region, and therefore provides one of the main constraints on those 
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processes. For a general review of the Kuiper belt size distribution, see lPetit et al.l 1)200. 
are three properties which describe the general shape of the Kuiper belt size distribution: 



There 



the existence of the largest objects, Eris and Pluto 

the large object size distribution for D > 100 km which is well chara c terized by a stee 
power-law ^ oc r~ q with q ~ 4.8 (jdadman et alJhoOll ; IPetit et alJbood : iFraser et alJboO 



the existence of a "roll-over" at r ~ 25 — 60 km where the size distribution flattens to a much 



shallower distribution tha n for larger objects ([Bernstein et al 
Fraser Kavelaarsl 12009) 



2004 



Fuentes &: Holman 



2008 



These three properties need to be reproduced by any successful attempt at recreating the accre- 
tionary and collisional history of the Kuiper belt region. While this has not yet been achieved, 
these properties provide some insight into the general history of the belt. 



Accretion simulations such as those of ISternl (|1996l ) ; iKenvon &: Bromlevl (|200ll ) and iKenvon 



(|2002l ) have demonstrated that objects as large as Eris could have accreted on timescales shorter 
than the age of the Solar system if the mass of the Kuiper belt was at least two orders of magnitude 
larger than the current mass, and if the relative encounter velocities, and hence, the eccentricities 
and inclinations of the objects were significantly lower in the early Solar system (see review by 



Kenyon et al.1 (j2008l )). The existence of the steep large object size distrib ution however, suggests 
that accretion was a short-lived process (jKenvonl 12002 : IFraser et al.ll2008l ). Some event must have 



disrupted accretion - likely the same event which scattered the primordial Kuiper belt onto orbits 
with high inclinations and eccentricities as observed today (see a review of proposed processes by 
Morbidelli et al.l (|2008l )) - otherwise the large object size distribution would be too shallow to be 
compatible with the observed distribution. 



The simulations of ISternl (jl996l ) ; IKenvon &: Bromley! (|200ll ) and IKenvon! (J2002J) cannot re- 
produce the large roll-over size. In these simulations, interaction velocities remain low, such that 
objects larger than r ~ 1 km do not experience disruptive collisions over the age of the Solar system. 
Thus, the size distribution for objects larger than r ~ 2 km exhibits a steep slope comparable to 
that observed today, and the accretion break, or roll-over, occurs at radii too small to be compatible 
with observations. This suggests that after the belt was dynamically excited, it remained massive 
enough for a time long enough to allow collisional erosion to produce the large roll-over observed 
today. 



A model presented by IKenvon Bromlevl (|2004l ) could reproduce the large roll-over if gravita- 
tional stirring by Neptune at its current location was included early in the simulations. Accretion 
occurred at a sufficient rate in these simulations to produce Eris-sized objects. In these simulations, 



1 Much of the observations that constrain the size distr ibution for radii r 
recent works I Fraser fc Kavelaars 20091 : Fuentes et al. 20091 ) . 



20 — 60 km are published in more 
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accretion lasted too long however, resulting in a large object size distribution too shallow to be com- 
patible with observations. In addition, the mass loss du e to collisional gr i nding was insufficient to 



produ ce the tenuous modern-day belt (see discussions bv lMorbidelli et al.l (|2008l ) and lKenvon et al 
(|2008l )). The Kuiper belt must have undergone a more rapid excitation and mass depletion than 
occurred in these calculations. 



Pan &: Saril (|2005l ) presented an anal ytical, order of m agnitude, collisional evolution model, 



reminiscent of the ground breaking work of iDohnanvil (| 19691 ) . With this model, they demonstrated 
that the large roll-over size could be produced by collisional grinding on timescales shorter than the 
age of the Solar system. They however, assumed that the size distribution rolled over to collisional 
equilibrium. It has been demonstrated that there is a range of objects with sizes smaller than 
the roll-over, which are pref erentially erod ed but which are not being replenished from fragments 
of larger disrupted objects (|Kenyonll2002l ). These objects do not achieve collisional equilibrium. 
Equilibrium is only achieved at a some sm aller radius, the ex act size of which, depends on the 
density of planetesimals. Thus, the model of lPan &; Saril (|2005l ) likely, does not predict the correct 
shape of the size distribution smaller than the roll-over, or the rate at which the roll-over evolves 
to larger sizes. 



Benavidez fe Campo Bagatinl (J2009J) present a collisional evolution model of the Kuiper belt, 
in which they include collisions from multiple dynamical classes. Using this model, they calculate 
the collisional evolution of the Kuiper belt, starting from an initially steep size distribution for 
al l sizes, with slope similar to that observed for large objects. Their results confirm the findings 
of lPan Saril (120051 ): collisional grinding with relative velocities comparable to that observed can 
disrupt the majority of objects with r ~ 50 — 100 km, and in effect, produce a roll-over at that 
size. These calculations however, likely do not reproduce the collisional history of the Kui per belt, 
as th ey start from an initial condition not expected from standard accretion processes (jKenvon 
20021 ). Accretion simulations which produce Eris-sized objects produce size distributions with an 
accretion break at r ~ 1 — 2 km, not a steep distribution for all sizes. 

Here we present a collisional evolution model, which we utilize to calculate the size distribution 
of the Kuiper belt from some starting condition. With this model, we wish to determine whether or 
not collisional erosion could produce the r ~ 20 — 50 km roll-over, starting from a size distribution 
similar to that produced by models of accretion in the outer Solar system, but in the dynamically 
excited conditions of the current day Kuiper belt. From these calculations, we wish to constrain 
the shape of the modern day size distribution smaller than the roll-over, and to constrain the mass 
of the Kuiper belt at the end of accretion. 

In Section [2] we present our collisional evolution model, and the planetesimal strength and 
shattering models we adopt. In Section [3] we present the main results of our calculation, namely 
that that existence of a break at r ~ 1 km left-over from accretion causes a divot in the size 
distribution at larger sizes. In Section 0] we present the consequences of our model, and finish with 
concluding remarks in Section [5j 
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The Model 



Our collisional evolution model uses a formalism similar to that of lOhtsuki et al.l (|1990l ) which 
uses bins of constant mass, but tracks average bin radiuJE The model considers a swarm of 
planetesimals distributed in a finite number of size bins, where bin i contains a number of obiects, 
Ni , with radii „ < r < „ + Sr., where S ~ 1.10 In our model we assuore that all particles 
have the same relative collision velocity, v re i. We do not include velocity evolution in our model, 
and assume that v re i is a constant during the simulations. Each collision can either result in 
catastrophic disruption and dispersal if the collision energy, Q, is sufficiently high, or cratering and 
mass accretion otherwise. We adopt the standard center-of-mass collision energy Q = }_Vih—i!h^L _ 
In the following sections we describe the details of the model and our prescription for the collisional 
outcomes. 



2.1. Disruption and Cratering Model 



We adopt the catastrophic disruption energy threshold functional form presented bv lBenz Asphaug 



(jl999j) . This form includes contributions from the tensile strength of the material body, as well as 



a gravitational binding term, and is given by 



Q*d = Qo(z — +Bp(- (1) 

V 1 cm/ V 1 cm/ 

where Q* D is defined as the disruption energy per unit mass of the target body required to shatter 
and disperse 50% of the target into a spectrum of fragments, p is the target density which we 
assume is constant for all sizes, and Q Q , a, b, and B are constants appropriate for the material 
properties of Kuiper belt objects. 



Benz &: Asphaugj (|1999l ) utilized smooth particle hydrodynamic simulations to determine the 
outcome of disruptive collisions between icy bodies. They calibrated Equation [1] under a large 
range of impact parameters (velocity, size ratio, etc.) for both basalt and water ice targets. They 
found that the disruption energy and hence the constants Q Q , a, b, and B depended on the relative 
impact velocity. They also found that the mass of the largest remaining fragment Mlrf depended 
linearly on the ratio of the impact energy to the disruption energy, Q/Q* D . They demonstrated 
that the mass of the largest remaining fragment can be well represented by 



2 Bins defined by an average radius rather than mass were adopted to ease analysis of the resultant size distributions. 
This choice however, has no effect on the results. 

3 Note that 5 here is the cube-root of the usual "mass-delta" - the ratio of masses in consecutive bins - in other 
calculations (see Kenvon fc Luu 19981 ) 
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Mlrf 
M t 



7 



Y-X 



Q 



M t Q* D 



(2) 



where Mt is the initial target mass. Examination of Figure 9. from lBenz fc Asphaus (|1999i ) reveals 
that Y = 1.0 and X = 0.55 in Equation [2] is an acceptab l e rep resentation of 7, and we adopt 
these parameters for ou r calculations. IStewart 8z Leinhardt] (|2009l ) have confirmed the findings of 
Benz &: Asphaua (|1999l ) and from their simula t ions have determined that X ~ 0.48. The small 
difference in X found by IStewart &: Leinhardtl (120091 ) and the value we adopt here will make an 
insignificant difference on our results, and will have no effect on our conclusions. 

We assume that the distribution of collisional fragments from a disruption is a power-law, 
and is given by dN/dr = Ar~ qD where qo is the logarithmic slope. It has been shown that in 
the asteroid belt, t he distribution of fragments is better represented by a two slope distribution 
(jBottke et alj|2005l ). In our calculations we considered a range of two sloped models. We found 
that, for a reasonable range of two-sloped models, the change in the size distribution results were 
small compared to the changes caused by adjustment of the strength and cratering laws - also much 
less than the variation introduced into the results from the unknown collisional history. Thus, we 
choose the power-law as a simple, scale free, one parameter representation that avoids including 
much of the uncertain physics of fragmentation in our model. 

Assuming that the density of planetesimals is constant for all sizes, from Equation [2] we see 
that the radius of the largest remaining fragment is given by 



r L RF = W 3 (3) 

where t% is the radius of the target body. Assuming there is only one largest remaining fragment, 
the normalization of the fragment model is given by A = (qp — V)^ qr> ~ 1 )/ 3 r^ to where qn > 1. 
Conserving mass, we have 

prLRF 

/ Ar~ qD m(r)dr = m(r^) (4) 
Jo 

where m(r) is the mass of an object with radius r. Equation U] can be solved for to give 

QD = ——■ (5 

7+1 

As 7 > 0, for our fragment model, 1 < qr> < 4. 

For impacts with energy insufficient to disrupt the target, we considering a cratering model 
where the excavated crater mass is given by 



Mcrat = OtffracQ 



(6) 
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where a is the crater excavation coefficient, which depends on the material properties of the target 
and impactor, and // mc is the fraction of kinetic energ y that went into shatt ering the target 
material. We consider values of a = 10 -8 — 10~ 9 s 2 cm -2 ( Petit fc Farinellalll993l ). 



We assume that the mass fraction of ejected material with velocity greater than v is given by 

-k 



fhi{v) 



V 



(7) 



where k = § (jPetit fc Farinellalll99.^ . By assuming that a fraction, Jke of the impact energy is 



imparted into kinetic energy of the excavated material (Equation [6]) , and requiring conservation of 

- ^ KE with the fragment mass escaping 



energy, the minimum ejecta velocity is given by v Tl 



9 Qtffrac 



the body given by M e . 



M 



rat 



where v, 



esa,j 



is the mutual escape 



velocity of colliding pair i, j (IWetherill &: Stewart! Il993l ). Clearly, accretion will occur if M esc is 
smaller than the mass of the impactor. 



We assume that the crater fragment distribution is a power-law with slope q c 
largest remaining fragment 



-3.4, and a 



rLRF,c 



4 + & 

47T/9 



M a 



rat 



1/3 



(8) 



Simulations of lBenz Asphaugi ()1999l ): lLeinhardt fc Stewart! ([2009J) and lStewart fc Leinhardt 

(|2009j) have demonstrated that catastrophic collisions still occur for energies as little as Q ~ 0.2 Q* d . 



If we consider a maximum crater size as a fraction of the target body radius, fcrat, an d require 
continuity between cratering and disruption regimes, the energy threshold for which collisions are 
just disruptive, as a fraction of Q* D , is given by fdis = y ~ ~ fcrat) 3 )- A minimum disruption 
energy criterion f^ 8 = 0.2 corresponds to fcrat = 0.04. In this way, events which strip more 
than ~ 10% of the mass of the target are considered disruptive collisions, while those that strip a 
smaller amount, are considered cratering events. In our calculations, we consider maximum crater 
sizes fcrat = 0.05 — 0.15 corresponding to minimum catastrophic disruption energy thresholds of 
0.4-0.7Q^. 

We note here that there are many effects other than size that determine disrupt ion strength and 



collisi onal outcome, such as impactor size and target porosity (see discussion by iLeinhardt et al 



(|2008l )). Given the uncertainty in the knowledge of the processes that shaped the Kuiper belt size 
distribution, the specifics of these effects are of little import, and we do not consider them in the 
collisional model given by equations [TKHJ 
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2.2. Collision Rate 

From simple cross-section arguments, it can be shown that, if objects from bin i are passing 
through the swarm of objects from bin j, the number of collisions between the objects from bins i 
and j, in time-step At is 



7r(r| + r])v rd At 



Ncoii^ = — ^ NiNj (9) 



where v re i is the relative encounter velocity, and V is the volume occupied by the swarm of particles. 

Because a single disruptive collision, or accretion of that object onto a larger target will remove 
an object from a bin, Equation [9] does not represent the true number of collisions between bins i 
and j - there can only be as many disruptive collisions as there are objects in a bin to be disrupted. 
Therefore, in calculating collision rates, one must consider disruptions of bin i by other impactors 
k / j, or if i is accreted onto a larger object. We adopt a probabilistic approach. We set the 
probability that an object from bin i is disrupted by an object from bin j (or accreted onto j if 
r j > r «)> Pi,jj as the probability that the disruption of the target i was from an impactor from bin 
j, PDij-, times the probability that an object from bin j was available to cause the disruption, Pj. 

Po i . is simply given by the number of collisions of objects from bin i by those from bin j, 
given by Equation [91 divided by the number of objects in bin i. That is Pd^ = — w^- 

The probability that impactor j is available can be found by setting Pj = 1 — Pj where Pj is 
the probability that an object j is disrupted by or accreted by an object with radius r& ^ fj. That 
is, Pj = 1 — -j^r Y2k^i Ncoll j k ■ The true number of disruptive collisions from bin i by bin j is then 

N disruptiJ =NiP Di .Pj. (10) 

Equation [TU] predicts a smaller number of collisions for large objects compared to Equation The 
difference imposses significant changes in the long-term collisional evolution when the number of 
disruptive collisions in a time-step is > 1% the number of objects in that bin. If the time-step is 
small enough, then the effect is small. 

For completeness we apply a corrective term to Equation [9l ( 1 + ( v ^"' j \ j to account for 



- Vrel 

gravitational focusing. For the impact velocities considered in this work, this term is negibibly 
small, and could be ignored without changing our conclusions. 
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2.3. Code Details and Nominal Parameters 

In any given time-step, the number of disruptive collisions per bin is first calculated with 
Equation [TUJ using an adaptive time-step, such that the number of objects being removed from a 
bin - either from disruption, or accretion and cratering - is no more than 1% the number of objects in 
the bin before the time-step. This minimizes code run-times while preserving the correct collisional 
evolution (see Appendix A). If the probability of collision from Equation [10] is less than unity, then 
a single collision occurrence is decided randomly, ie. the collision occurs if R < NiPuijPj, where 
R is a randomly generated number with < R < 1. 

Formally, if the number of collisions of a bin is less than 10 8 per time-step, the number of 
collisions is kept as an integer va lue, and a random n u mber generator is used to determine if the 



number is rounded up or down. IWetherill &: Stewart] f| 19891 ) has demonstrated that requiring an 



integer number of collisions per time-step when the number is less than 10 8 is a sufficient condition 
to produce the correct time averaged collision rates. 

For each non-disruptive collision between a large object i and a small object j, the resultant 
body has a mass m(rj) + m(rj) — M esc . An object is removed from bin j, and if the resulting mass 
is sufficiently large, a body is removed from bin i and added to the bin with the appropriate bin 
k. For these collisions, their is a difference in mass between the resultant body, and the average 
mass of the bin it was promoted to. This mass, accreted by bin k, M acc ^ (negative if cratering 
occurs) is accumulated over multiple time-steps. When the accreted (or cratered) mass becomes 
large enough, objects can "grow" (or shrink) to subsequent bins, ie. when M acc ^ > Mjy-i — 
(sign set according to accretion or cratering) objects are moved randomly from bin % to k ± 1. This 
is done in such a way as to reflect the fact that the accreted mass is occurring over all objects in 
bin k and not just one. Thus, \ M^±°-M k \ objects are added to bin k± 1 from bin k when M ^ c,fc > R 
where R is a randomly generated value with < R < 1. This process preserves total system mass. 



We consider bins increasing logarithmically in radius. That is rj = rj_i<5 where 5 > 1. For our 
standard simu lations, we con sider 5 = 1.1. It has been shown that a finite 5 produces an acceleration 



in the results (jWetherilll (see I199Q ). and Appendix A). THus, our conclusions drawn here have the 



caveat that the evolution might occur in nature slightly slower than in our simulations. 



Campo Bagatin et al.l ([1994) has demonstrated that a finite minimum size bin will cause sub- 
stantial wave-like deviations of the small size distribution away from the expected collisional equi- 
librium distribution; the smallest particles have no smaller impactors, and are in over abundance 
compared to larger particles. This over abundance creates an increased collision rate for larger 
particles, and so on. To prevent such a behavior, we enforce a collisional equilibrium slope for 
the smallest bins by removing excess objects from those bin s. While this meth od is not the most 



accurate way of correcting for the finite bin minimum (see iKenyon et al.1 120081 . for a discussion) , 
it is trivial to implement, and ensures accurate enough small object treatment to not affect our 
conclusions about the > 1 km size distribution. It was found that enforcing collisional equilibrium 
over the smallest 30 — 40 bins is necessary to prevent the over abundance of the smallest objects 
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(see Appendix |A|) . 



Our calculations were performed on an NVIDIA GTX 280 graphics card using the CUDA 
programming environment Q These cards provide 30 processors, each containing 8 processor cores, 
all running in parallel, and excel at highly parallel calculations, with speeds up to 50 times higher 
than typical multi-core CPUs. Due to the nature of the calculations however, we are limited to 
2 n bins. For our nominal calculations, we use 256 bins with 5 = 1.1 for our nominal simulations, 
allowing us to probe more than 10 orders of magnitude in radius. We considered larger 5 with 
fewer bins to determine the effects of bin width and minimum bin size on our results. We found 
that changing the bin size did not produce a noticeable affect on our conclusions considering the 
range of results caused by the uncertain shattering and cratering models. Thus, we adopt 5 = 1.1 
for all the calculations we present here. 

We wish to model collisional evolution that occurs in dynamical conditions typical of the 
modern day Kuiper belt. To that extent, we roughly model the Kuiper belt as a single annulus, 
with inner and outer radii of 30 and 60 AU. We set a scale height of 20 AU similar to the scale 
height of the current Kuiper belt. We note th at this scale he ight is significantly larger than that 
achieved in standard accretion calcuations (see iKenvonl 120021 . for discussion) . Such a scale height 
can only be achieved from scattering dynamics not included in those acc retion models, and wh ich 
is necessary to shape the primordial Kuiper belt to its current state (see iMorbidelli et al.l 1200a ). 



It has been shown that accretion calculations which util ize a single a nnulus, as done here, 
produce different results than those which utilize many annuli (|Kenvonll2002l ). The calculations we 
present here are intended to show the collisional evolution of an excited Kuiper belt. Standard multi- 
annulus codes cannot accurately represent the dynamical struc ture of the Kuiper belt and therefor e 
do not accurately model the collisional evolution of this region (jBenavidez &: Campo Bagat in 2009). 
Rather, a proper treatment w ould involve a combined N-body a nd collisional evolution calculation, 
as done in lBottke et al.l (|2005l ) and lCharnoz Morbidellil (|2007l ). Given that the dynamical history 
of the Kuiper belt is not yet known, such a calculation is beyond the scope of this work. As such, 
the calculations presented should only be interpreted as rough approximations to the true collisional 
evolution that occurred in the Kuiper belt. 



We set v re i = 1 km s _1 which iDell'Oro et al.l ( 20011 ) has demonstrated is a typical collision 
velocity in the modern Kuiper belt. At this velocity, targets of r > 5 km will be disrupted in 
collisions with size ratio larger than ~ 1 : 10 for the impactor and target. Collisions with smaller 
impactors will typically cause cratering, and accretion will occur - slowly - for only the largest 
(r > 250 km) objects. 

For our simulations we define two strength models. For the "strong" mod el, we adopt the 
streng th parameters of water-ice for a 0.5 km s -1 collision velocity presented by Benz fc Asphaug 



dl999j), withQ = 7xl0 7 erg g _1 , B = 2.1 erg cm 3 g" 



0.45, and b = 1.19. Leinhardt Sz Stewart 



4 www.nvidia.com / cuda 
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(|2009l ) however, has demonstrated that the above destruction criterion is likely a factor of ~ 3 
too high for objects i n the Kuiper belt. Thus, we adopt the strength parameters suggested by 
Leinhardt &: Stewart! (|2009l ). which have the same slopes a and b, as the strong model, but with a 
factor of 3 less disruption energy at all sizes for our "weak" model. We utilize the strong model for 
our nominal simulations, and discuss the effects of using the weak model in Sections [3] and 01 



= 0.05, and part ition the collision e nergy into fracturing 
0.1. Following iPetit &: Farinellal (|1993l ). for our nomi- 



We set the maximum crater size, fcrat 
and kinetic energy as ff rac = 0.8 and Jke : 
nal model, we set the cratering velocity fragment k = |, the cratering efficiency a = 10 -9 s 2 cm -2 
corresponding to "strong" surfaces or low cratering efficiency, and the fragment size distribution 
slope q c = 3.4. 



2.4. Initial Conditions 



The goal of this work is to determine the collisional evolution of the Kuiper belt size distribution 
after the epoch o f accretion. To that end, we consider a "typical" size distribution found from accre- 
tion simulations (jKenvon Luulll998l ; iKenvon &: Bromlevll200ll . 12004 ; ISternlll996l ; IStern &: Colwell 
19971 ). That is, one which has a steep power-law for the largest o b jects with qi > 4. At some smaller 



the slope becomes shallower, o r breaks to q2 ~ — 1 to 1 (jStern 



199 



Stern fc Colwell 



1997 



Kenvon|[200^ : IKenvon Bromlevl[20o3 ) - we utilize a sh arp break for simplicity The "break-radius" 
r bi - typically 2 km for the Kuiper belt (|Kenvonll2002l ). corresponds to the size that objects which 
are smaller have undergone at least one disruptive collision. At some even smaller size, r&2 ~ 0.5 km, 
collisional equilibrium is reached, and the size-distribution slope turns up to the canonical q% = 3.5. 
Again we assume a sharp slope-change. The collisional equilibrium slope ty pically only deviates a t 
sizes small enough such that radiation effects remove objects from the belt (IThebault et al.l 12003). 
The effect however, is inconsequential on our results, so we do not consider radiation effects here. 

The most recent measurements of the Kuiper belt size distribution have determined that q\ ~ 



4.8 (Fuentes &: Holman 



2008 



Fraser &; Kavelaarsl l2009l) and breaks at a radius 



The steep slope implies a short accretion timescale (Kenyon 



been preserved from the epoch of accretion (jDurda &: Stern 



2002 



200C 



Fraser 



et al 



Pan &: Sari 



20 - 60 km. 



20081 ) and has 



20051 ). Thus, for our 



nominal simulations, we start with an initial large object slope equal to that observed in the modern 
day b elt, ie . qi = 4.8. We start with break radii and slopes inspired by the simulations of [Kenyon 
(J2002J) and Kenyon fc Bromley! (J200J) which represent the most complete simulations of accretion 
in the Kuiper belt range to-date, and for our nominal simulations we take qi = 0.0 — 2.0, q^ = 3.5, 
r^i = 2 km, and r&2 = 0.5 km. We discuss variations of these initial conditions in Section UJ 

Accretion simulations of the region ha ye demonstrat ed that the Kuiper belt must have had 
a few orders of magnitude more mass (s e e iKenvonl 120021. for a review) than currently o bserved 
(|Trujillo et al.l l200ll : Idadman et al.l I2OO1I : bernstein et al.l I2OO4J : iFuentes Holmanl hood ) , other- 
wise Pluto could not have achieved its size over the age of the Solar system. The exact initial mass 
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available at the late stages of accretion however, is unknown. For our nominal simulations we adopt 
a 45 M ffi Kuiper belt and vary this mass to test the evolution for different planetesimal densities. 
We note her e that 45 Mm i s at the large end of the mass thought to have existed in the primordial 
Kuiper belt (|Kenvonll2002l ; iGomesI 120031 ) . This, the rate of evolution produced in calculation with 
such a high mass should be interpreted as upper limits to the real evolution. The parameters for 
our nominal simulations are listed in Table [TJ 



3. Results 

In Figure [TJ we present the results of our nominal simulations at specific times. The expected 
behavior of the collisional erosion occurs at the large and small scales. For particles smaller than 
the initial r^, because these particles were initialized with the collisional equilibrium slope, the 
overall trend is to remain with that slope, and to decrease in mass as more particles are ground 
away. For the largest (r > 150 km) objects, very few (if any) disruptive collisions occur. Accretion 
and cratering results in changes in the large object size distribution by less than ± ~ 1% over the 
age of the Solar system. 

For objects with r ~ 10 a distinct divot forms in the size distribution. This behavior is 
caused by the initial break or turn-over at r^i) the decrease in slope causes an absence in the number 
of disruptors available to shatter objects with r ~ r^x- Thus, due to the change in slope at r^x, any 
objects capable of being shattered by impactors with r ~ r^i experience an enhanced collisional 
disruption rate. This in turn reduces the disruption rate of objects capable of being shattered by 
the impactors with r ~ 10 r^i, which produces a distinct kink where the size distribution rolls away 
from the initial accretion slope q\ into the divot. 

The location of the divot can be roughly estimated by equating the collision energy of impactors 
at the accretion break radius to the disruption energy of the bodies with the divot radius r^. 
Ignoring the tensile strength term in Equation (TJ and solving for r^iv we find 

r div = (Af dis Bp)W> v£s r^. (11) 

Equation [TT1 demonstrates that given the accretion break size r&i, the radius of the divot, r d i v , 
depends on the strength parameters of the gravity term in Equation [TJ B, and b, the density of 

Kuiper belt objects, p, the average interaction velocity, and the energy at which impacts transition 

-i 

from the cratering regime to the disruption regime, given by f d i s - The factor {Af d i s Bp) 3 + b however, 
is of order unity. Thus, the central divot radius depends mainly on the slope of the gravity strength 
scaling-law b, and the velocity at which objects interact. These conclusions are confirmed in Fig- 
ured! where we present the results of our nominal simulations with different strength parameters. 

The rate at which the divot deepens depends on the slope change at r&i; the greater the 
difference between q\ and q2, the greater the rate at which the divot deepens (see Figure [p. The 
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divot formation timescale can be roughly estimated from Equation [9j The divot is formed primarily 
from catastrophic collisions of objects larger than the initial break radius r&i where the initial size 
distribution is given by a power-law with slope q%, ie. ^ = Cr~ Ql where C is a normalization 
constant. Thus, by inserting a power-law for the number of smaller objects Nj in Equation [9l 
integrating over all sizes of object capable of disrupting objects at the divot, ie. from r&i to r^, 
and solving for time, we find that the timescale for divot formation t^ v is given by 



V 



N Tdiv J nvCr° div 



3-<?i 



1 _ f l ~1^ 1 — Z 3- "? 1 



l-qi 



3 - qi 



(12) 



where 



JV rj . 



is the fraction of objects with radius r^ v that have been disrupted, and / = ~ 0.2 



(see Equation [TTj) . By adopting the nominal parameters, and saying a divot has formed when yyyo 
of objects at the divot have been disrupted, ie. ^ c = 0.99, then the divot formation timescale 

r div 

is tdiv ~ 25 Myr. For this timescale estimate, we have ignored disruptions of objects with radii 
r bi < r < r div Thus, Equation [12] represents a lower-limit of the true divot formation timescale. 
The true divot formation timescale is ~ 50% longer if the break slope is qi = and ~ 3 longer if 
q2 = 2 (see Figure [1]) . 

Since the probability of a single object having a collision is approximately linearly dependent on 
the density of particles, the rate at which the divot forms is also linearly dependent on the density of 
particles (see Equations [91 and I12[) . Minimal variation in the results occur due to order of magnitude 
variations in the belt density at equivalent points in the simulations. This is demonstrated in 
Figure El where we compare the nominal simulation with a simulation using the same parameters 
with the density dropped by a factor of 10. 

Presented in Figure 0] are the results of our nominal simulations when the parameters relevant 
to the cratering model are changed. With certain strength models, namely those where cratering 
is responsible for removing a large fraction of object with r > lOr^i, small waves can form at 
those sizes. But for the majority of the simulations, the roll-over is smooth up to the edge of 
the divot where the density of objects decreases very rapidly. Other than the maximum cratering 
radius parameter f^ s which affects the divot location, little variation in the size distribution for 
r > rt,2 results from variations in the cratering parameters. Changes do occur for r < r&2- As there 
currently exists no measurement of the shape of the Kuiper belt size distribution in this size range, 
we leave the study of these parameters to later works. 



4. Discussion 



The divot produced in our calcul ations could provide a natural explanati on for the observed 
roll-over at r ~ 25 — 60 km dete c ted bylBernstein et al.l (j2004l ) and confirmed bv lFuentes fc Holman 



(|2008l ) and lFraser fc Kavelaarsl (|2009l ). That is, the size-distribution deviates away from the steep 
large-object accretion size-distribution into a divot at r^iv ~ 10 — 15 km. If the observed roll-over 
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is caused by a break left-over from accretion, then, from Equation [TT\ we find that the accretion 
break must have occurred at a radius, r&i ~ 0.5 — 3.5 km, for our adopted strength scaling slope 
b = 1.19 . The implied accretion size break r^i is consistent with that found in previous accretion 
models (jKenvonl 12002 ). 



Comparison of the ratio of objects larger and smaller than the kink radius, rj~, from the 
simulations, to the roll-over observed in the Kuiper belt can provide constraint on the mass of the 
belt at the end of accretion. Assuming that the accretion break-slope is, q2 ~ 0, and that the 
Kuiper belt strength scaling is well represented by our assumed value b = 1.19, a 45 M ffi Kuiper 
belt would produce a divot consistent with observations in > 40 Myr if r^i = 2 km and > 100 Myr 
for r&i = 1 km. r^i =2. The divot formation timescales are a fa ctor of ~ 2 longer for the in itial 
accretion break slope q<i = 2 and ~ 2 shorter for the scaling law o f Lcinh ardt h, Stewart! (|2009i ). ie. 
if objects are a factor of 3 weaker than suggested by Benz Sz AsphaugT f 1999 ). The time estimates 
also scale linearly with the density of planetesimals at the epoch of the dynamical excitation. 

Our Kuiper belt mass estimates are consistent with accretion simulations which suggest that 
the primordia l Kuiper belt had a mass > 30 M e (jStern Colwelll 119971 ; iKenvon Luul Il998l ; 
Kenvonl |2002j) . If we assume that the mass depletion could not have occurred more than 1 Gyr 
after the Kuiper belt was excited, then the minimum Kuiper belt mass required to produce a divot 
compatible with the observations is ~ 1 — 5 M®. 



The simulations of lCharnoz &: Morbidellil (120071 ) are the first attempt at coupling the dynamical 
and collisional history of the Kuiper belt. Using a similar initial mass distribution (45 M® between 
5 and 50 AU) as we have adopted in our simulations, they found that the mass depletion which has 
produced the low-density modern day Kuiper belt must have been a result of combined collisional 
grinding and dynamical effects. Otherwise , the planetesimal po pulations external to the Kuiper 
belt - the Oort cloud and scattered disk (see lGladman et al.ll2008l ) - would be too anemic by at least 
an order of magnitude. They found that the mass depletion of the Kuiper belt occurred over ~ 100 
Myr. They however adopted an initial size distribution that broke from a steep large-object slope 
to that of collisional equilibrium rather than the accretion model inspired distribution we adopt 
here, complicating comparison of the two works. Roughly speaking however, the mass depletion 
timescale they found would be sufficient to produce a divot, if starting from an accretion-like size 
distribution. They also found that populations external to the Kuiper belt underwent significantly 
enhanced collisional evolution compared to the Kuiper belt. This suggests that the scattered disk 
and Oort cloud populations would have very few objects with 1 < r < 20 km. 

Our results are consistent with some of the peculi arities of detected K uiper belt objects. Only 
one object with radius r < 15 km has been detected (jFuentes et al.l l2009). Such an observation is 
consistent with a flat q ~ power-law for objects smaller than the current roll-over. It is suggestive 
however, that there is a substantial depletion of objects with similar size. This is consistent with 
the idea that the size distribution is not a power-law for r ~ 10 — 20 km, but rather has a shape 
like the divot produced in our calculations. The shape of the luminosity function presented by 
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Fuentes et al.l (j2009l ) hints that a divot is detected at m(R) ~ 26.5. The observational evidence of 
this however, is extremely sparse. More observations of small r ~ 10 — 40 km Kuiper belt objects 
are required before the true size distribution shape is known. 

We find that objects with sizes similar to the initial accretion break r^i do not suffer signifi- 
cant depletion on timescales necessary to form the divot. Our findings demonstrate that the size 
distribution should exhibit a depletio n of objects for r ~ 10 km compared to objects with r ~ r^. 



Comparison of our results to that of iBenavidez &; Campo Bagatinl (|2009l ) demonstrates that the 



shape of the size distribution caused by collisional evolution is highly dependent on the initial size 
distribution; the break can be formed with different initial conditions, but the location, and inter- 
pretation of the break change. Our results demonstrate that the detection of a divot in the size 
distribution is distinct evidence for the existence of a break in the size distribution at the onset of 
collisional disruption. 

Our calculations demonstrate that, with a detection of the remaining peak at r&i, and the divot 
center r^iv, the strength scaling slope b as well as the accretion break radius r^i can be inferred. This 
would provide a stringent test of numerical impact simulations of icy-bodies, as well as a strong 
constraint on accretion models of the Kuiper belt region. The current uncertainty on the roll- 
over size cannot place anymore constraint on accretion or shattering physics of icy planetesimals, 
other than to say that the observations are consistent with the current understanding of these two 
processes. 

Our findings suggest that the sky density of objects with r ~ 1 km is ~ 10 6 - 10 7 objects per 
square degree. Detection by serendipitous occultation of star light is currently the only method 
of detection sensitive to these small KBOs. These exciting exp eriments have curren t ly placed an 



uppe r limit of ~ 10 8 — 10 9 objects per square degree at this size (jBickerton et al.ll2008l ; IZhang et al 



2008) and future experiments will ultimately determine if the roll-over is caused by the mechanism 



we have suggested here. 

Our calculations, and the existence of the steep slope observed for the large object size dis- 
tribution, suggests that the Kuiper belt underwent a short period of accretion, where objects as 
large as Eris were produced. Some dynamical event truncated accretion by increasing interaction 
velocities to disruption. On a short ~ 10 — 100 Myr timescale, a divot, or depletion of objects 
with r ~ 10 km was produced before the dynamical scattering removed most of the mass from the 
system, freezing the shape of the size distribution, which has undergone minimal evolution to the 
current day. 



5. Conclusions 

Our calculations have provided the first direct link from the "typical" size distribution produced 
from accretion, to the collisionally modified size distribution of the modern day Kuiper belt. Our 
calculations have demonstrated that, given a size distribution similar to that produced by accretion, 
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the likely result of collisional evolution in the Kuiper belt, is to produce a divot in the Kuiper belt 
size distribution at r ~ 10—15 km, and a roll-over from the large object accretion slope at r ~ 25—60 
km, compatible with the observed roll-over in the Kuiper belt size distribution. We conclude that 
the size distribution is likely not a power-law for objects smaller than the roll-over, and exhibits a 
dearth of objects at the divot location. Assuming a divot has been formed, then a measurement of 
the divot center radius would provide constraint on the strength scaling law in the gravity regime 
for KBOs, and constrain the size distribution of r ~ 1 km objects left-over from the accretion 
processes in the early Solar system. 

From our calculations, we find that the Kuiper belt must have had at least 1 — 5 M ffi of material 
after its objects were excited onto orbits with eccentricities comparable to those observed, otherwise 
the roll-over would not be produced on Gyr timescales. If the belt had a mass of ~ 45 M©, then a 
divot would form in > 40 Myr. 

From the results of our calculations, we predict a deep absence of objects at r ~ 10 km. We 
also predict that the sky density of objects with r ~ 1 km is ~ 10 6 — 10 7 per square degree, 2-3 
orders of magnitude larger than if the size distribution is a s hallow power-law for objects smaller 



than the roll-over, as suggested bv lFraser k, Kavelaard (|2009i ). 
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A. Appendix A 



Here we present va rious tests to confirm the correct behavior of the collisional model. As done 
by Qhtsuki et al. ( 1990 ) and Wetherii] ( 199C ). the behavior of a collisional evolution model -at least 
in terms of accretion - can be assessed by comparing numerical simulations of coagulation to the 
analytic solutions of the discrete coagulation equation given by 
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^t Uk = \ ^ Ai 3 nin 3 ~ U k^2 A ik n i- ( Al ) 

i+j=k i=l 

Here, j^n^ is the time rate of change of objects in bin k from collisions adding objects to (left 
term) and removing objects from (right term) bin k summed over all bins. A^ is the "kernel", and 
governs the probability of collisions between i and j. This is usually written in terms of mass bins % 
and j rather than the size bins we consider above. Any successful accretion or collision model must 
reproduce analytic solutions of Equation IA11 Discrepancies between the numerica l and analytic 



soluti ons can highlight possible caveats about the results of the numerical models. lOhtsuki et al 



(|1990l ) has demonstrated that, in his constant mass bin formalism, the numerical solutions are 



typically accelerate d compared to th e analytic counter parts. It is worth noting that, in the mass- 



batch formalism of IWetherilll (|1990l ). the numerical solution lags the analytic one. As we use the 
constant mass bin formalism, the conclusions about the behaviour of his model are applicable here. 
We discuss this below. The acceleration decreases with the bin-spacing 50 and varies with each 
kernel used. Thus, comparison of the numerical and analytic solution can provide a metric for the 
largest allowable 5. These comparisons also provide a measure of how large a time-step can be used 
before the incorrect evolution is produced. 

The simplest kernel is a constant, ie. = A = constant, ie. the collision probability doesn't 
depend on target mass. This kernel produces smooth orderly growth, and the solution to the 
coagulation equation for this kernel, ie. the number of objects in a bin, njt, is ^ = f 2 (1 — f ) k l 
where / = (1 + r//2) _1 and n = \n D t is the dimensionless time. n Q is the total number of objects 
at time zero, all of which are of the smallest mass. 

The analytic solution and the numerical solution of the coagulation equation with the constant 
kernel are presented in Figure [5] utilizing a 5 = 1.1, n a = 10 18 particleqj, 128 bins, and by not 
allowing the number of objects in a bin to change by more than 0.1% per time-step. As can 
be seen, the numerical s olution matches the analytic solution well. The discrepancies present are 



similar to those found by lOhtsuki et al.1 (J1990). The numerical solution produces more large objects 



than in the analytic one. The difference however, is almost negligible and increases for larger 5. 

Another kernel for which an analytic solution exists is for collision probabilities proportional 
to the sum of colliding masses, ie. Aij = (5{i + j). The analytic solution for this kernel is = 
n Q k k t f (1 — f) k ~ l e~ k ^~^ where / = e~ v and n = (3n t is the dimensionless time. The analytic 
and numerical solutions for this kernel are presented in Figures [6] and [7] for n a = 10 18 particles, 
128 bins, and 5 = 1.080 This kernel creates large objects more rapidly than the constant kernel. 



5 As done in the main text, we still refer to S here in terms of object size. Our 5 is the cube of that used in the 
coagulation equation solution. 

6 it was found that for n > 10 10 , the results were independent of the choice in n 

7 chosen to ensure the largest mass bin is larger than the largest object in the analytic solution. 
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As such, the numerical solution is more sensitive to the choice in 5 than for the constant kernel. 



As shown by lOhtsuki et al.l (|1990l ) , a finite 5 in the numerical solution creates an acceleration in 
the coagulation, as well as an over-abundance of large objects compared to the analytic solution. 
Our model exhibits virtually identical behavior. We find that the acceleration is nearly a constant 
over the simulations (see Figure [7|). That is, for 8 = 1.08 the shape of the numerical solution best 
matches the analytic solution at a time ~ 12% earlier than expected. This factor is constant over 
the simulations, and scales with the choice in 5. 

The kernel proportional to the product of the colliding masses, ie. Aij = 7 (i * j) produces 
runaway growth at rj = 1 (again rj = jn a t is the dimensionless time). As such, the kernel is difficult 
to solve analytically, as the time-step must vary with mass of the run-away body to re-create the 
proper evolution. The analytic solution for the non-runaway bodies is n/% = " o( -^ — 1 e~ kv , 

and the mass of the runaway body is m = n exp 

\-JJ2k=i k 2 ^dr]'} (|Wetherilllll99(ih . 



The code presented here was intended for collisional evolution calculations, where the primary 
process is disruption of objects to ever smaller sizes. The code inefficiently handles the case of 
runaway growth; the code takes many more steps than necessary, causing the runtime to increase 
substantially as runaway continues. Our code was run only long enough to demonstrate that the 
numerical solution produces the correct mass of the runaway body, and distribution of smaller 
bodies (see Figured]). Future versions of the software will include modifications to handle runaway 
growth more efficiently. But as the simulations presented above have minimal accretion of even the 
largest objects, inclusion of these modifications was unnecessary, and omitted for this work. 

Presented in Figure [8] is a comparison of the numerical and analytic solutions to the coagulation 
equation before and after runaway begins using n a = 10 18 and 5 = 1.06, and allowing the number 
of objects per bin to vary by 0.1% per time-step. Much like the previous examples, the numerical 
solution produces an accelerated evolution compared to the analytic solution. The acceleration in 
the numerical solution of runaway growth is not constant with time, as is found in other numerical 



models (see for instance iKenvon &: Luul ll998). The acceleration is most pronounced for short times 



(~ 20% when rj = 0.4), and decreases with time. 

In the numerical solution, runaway growth occurs at rj = 0.976, slightly early compared to 
the analytic solution. The post-runaway numerical solutions are presented in Figure [8] at rj = 0.98 
and 1.01, and very accurately match the analytic solutions at rj = 1.0003 and 1.0006 respectively. 
The mass of the runaway body in the numerical solution at these two times, is 3.5 x 10~ 4 and 
5.6 x 10~ 4 the total system mass. The mass of the runaway body in the analytic solution at the 
times where the analytic and numerical solutions match is 3.5 x 10 -4 and 6 x 10 -4 the total system 
mass, almost identical to the numerical solution. The results presented in Figure [8] demonstrate 
that our code accurately reproduces the analytic solution to the coagulation equation for the case 
of runaway growth, before and after runaway commences. 

The solution to this kernel is much more sensitive to the choice in 5. For 5 > 1.1, the 
acceleration verges on unacceptably large, and the mass of the run-away body is not correctly 
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reproduced. Therefore, 5 = 1.1 is an approximate upper-limit on the bin spacing for simulations in 
which runaway growth occurs, and which utilize our model. 



B. Appendix B 

Here we present tests of the collisional evolution considering various details of our model. 
Specifically we test the effects of time-step on our collisional evolution calculations and the numerical 
solutions to the coagulation equation. Additionally, we test our adopted collisional equilibrium 
forcing over the smallest bins. 

The choice in time-step modifies the results of the calculations; too large and the results will not 
accurately produce the correct evolution. While large time-steps may repeat the general behavior 
of simulations which use smaller time-steps, significant stochastic variations in the results occur if 
the time-step used is too large. It was found that not allowing the number of objects in a bin to 
change by more than 1% per time-step was sufficient to ensure accuracy in our collisional evolution 
simulations as well as the solutions to the coagulation equation for kernels A = A and A = (3{i + j) 
(see Figures El El and [9]) . Further tightening this condition did not change the results by more than 
a few percent, while relaxing this condition further produced incorrect results for all numerical 
solutions. Note: none of the numerical solutions of the coagulation equation include a probabilistic 
treatment of the collision rates as given in Equation [10J Therefore, these solutions suffer the most 
from larger time-steps than do the main simulations presented here. 

For the simulations in which collisional equilibrium was forced for the smallest bins, we tested 
over how many bins this should be applied. The results are shown in Figure [TUl where we present 
the evolution of a system already in collisional equilibrium, but forcing equilibrium on the smallest 
30, 40, and 50 bins. 

If equilibrium is forced over 30 bins only, a large mass build-up occurs at the small-size end. 
This build-up increases in magnitude over time, and produces unrealistic results. If collisional 
equilibrium was forced over 40 and 50 bins, the mass build-up no longer occurs, and the smallest 
object become collisionally depleted, as expected. The minimum number of bins required to prevent 
the mass build up varies with the cratering and fragment distributions chosen for the simulations. 
The rate of depletion does depend on the number of bins over which collisional equilibrium is forced. 
The results however, are similar for objects more than ~ 2 orders of magnitude larger than the 
smallest bin. This effect will not affect the main conclusions drawn above about the large object 
size distribution. 



As demonstrated by lThebault et al.1 (j2003l ). proper treatment of the smallest objects requires 
modeling of the radiation forces which are ultimately responsible for removing the dust out of 
the region. Such a treatment is beyond the scope of this work. We choose to force collisional 
equilibrium over the smallest 40 bins, as a compromise between having an unrealistic small object 
mass buildup, and potentially removing mass too quickly from the system. All conclusions drawn 
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above have the caveat, that the small object size distribution might not be correctly accounted for 
here. 
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Table 1. Nominal Simulation Parameters 



Parameter 


Adopted Value 


Disruption Parameters 


Qo 


7 x 10 7 erg g 1 


B 


2.1 erg cm 3 g~ 2 


a 


-0.45 


b 


1.19 


X 


0.55 


Y 


1 


Craterinj 


I Parameters 


a 


10~ 9 s 2 cm -2 


Ike 


0.1 


ffrac 


0.8 


fCrat 


0.05 


<?c 


3.4 


k 


9 
4 


Miscellaneous Parameters 


P 


1.2 g cm 3 


Vrel 


1 km s _1 




0.01 cm 


5 


1.1 


Initial Size Distribution Parameters 


qi 


4.8 


Q2 


0-2 


Q3 


3.5 


m 


2 km 


H2 


0.5 km 
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Fig. 1. — Results of the nominal simulations with the accretion break, r&i = 2 km at different times 
in the evolution. Dotted line - initial distribution. Dashed line - 25 Myr. Dashed-dotted 
line - 100 Myr. Solid line - 500 Myr. a q<i = 0. b qi = 1. c q<i = 2. The arrows mark the 
approximate divot center radius r^ v predicted from Equation [TTJ NOTE: The radius bin-widths 
increase proportionally with the bin-centers. This correctly presents the number of objects in a 
bin, but causes the size distribution slopes to appear —1 shallower than the true slope. 
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Fig. 2. — Results of the simulations with the accretion break, = 2 km at 500 Myr but varying 
parameters relevant to the location of the divot radius rdiv The arrows mark the approximate 
divot radii predicted from Equation [Til 
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Fig. 3. — Results of the simulations with the accretion break, r&i = 2 km at 50 Myr but varying the 
density. Top: Solid line - nominal parameters with 45 M ffi , Dashed line - nominal parameters 
with 4.5 M®. Bottom: ratio of the two simulations presented above. 
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Fig. 4. — Results of the simulations with the accretion break, r^i = 2 km at 500 Myr but varying 
parameters relevant to cratering. 
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Fig. 5. — Analytic (solid line) and numerical (points) solutions to the coagulation equation at 
different time intervals for the kernel A = A. 
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Fig. 6. — Analytic (solid line) and numerical (points) solutions to the coagulation equation at 
different time intervals for the kernel A = (3(i + j). Dots and squares represent results when the 
bins are allowed to vary by 0.1% and 1% respectively. 
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Fig. 7. — As in Figured! but with the numerical solution presented at a time 12% earlier than the 
analytic solution. 
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Fig. 8. — Analytic (solid lines) and numerical (points) solutions to the coagulation equation 
at different time intervals for the kernel A = 7(27). Left: Before runaway growth commences. 
Presented are numerical solutions at 77 = 0.4, 0.7, and 0.976, along with analytic solutions at 
77 = 0.4, 0.7, 0.976, and 1. Right: Snap shot of the large non-runaway bodies, where the evolution 
is most rapid, after runaway growth commences. Presented are numerical solutions at r\ = 0.98 
and 1.01 along with analytic solutions at rj = 1.0003 and 1.0006. 
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Fig. 9. — Top: Results of collisional evolution for the nominal parameters at 100 and 500 Myr when 
allowing the number of objects in a bin to change by 0.1% (solid) and 1% (dotted) per timestep. 
Bottom: Ratio of the simulations presented in the top using two separate timesteps at 100 (dashed) 
and 500 (solid) Myr. 
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Fig. 10. — Results of collisional evolution after 2000 time-steps, starting from a distribution in 
collisional equilibrium (dotted line). The various lines are the results when equilibrium was forced 
over 30, 40, and 50 bins. 



